#include "base_funcs.h"

void display(gsl_vector* x) {
	for (int i = 0; i < x->size; i++) {
		cout << gsl_vector_get(x, i) << " ";
	}
	cout << "\n" << endl;
}

gsl_vector* add(gsl_vector* a, gsl_vector* b)
{
	gsl_vector* v = gsl_vector_alloc(a->size);
	gsl_vector_memcpy(v, a);
	gsl_vector_add(v, b); // v = a + b
	return v;
}

gsl_vector* add(gsl_vector* a, double b)
{
	gsl_vector* v = gsl_vector_alloc(a->size);
	gsl_vector_memcpy(v, a);
	gsl_vector_add_constant(v, b); // v = a + b
	return v;
}

gsl_vector* add(double a, gsl_vector* b)
{
	gsl_vector* v = gsl_vector_alloc(b->size);
	gsl_vector_memcpy(v, b);
	gsl_vector_add_constant(v, a); // v = a + b
	return v;
}

gsl_vector* sub(gsl_vector* a, gsl_vector* b)
{
	gsl_vector* v = gsl_vector_alloc(a->size);
	gsl_vector_memcpy(v, a);
	gsl_vector_sub(v, b); // v = a - b
	return v;
}

gsl_vector* sub(gsl_vector* a, double b)
{
	gsl_vector* v = gsl_vector_alloc(a->size);
	gsl_vector_memcpy(v, a);
	gsl_vector_add_constant(v, -b); // v = a - b
	return v;
}

gsl_vector* sub(double a, gsl_vector* b)
{
	gsl_vector* v = gsl_vector_alloc(b->size);
	gsl_vector_memcpy(v, b);
	gsl_vector_add_constant(v, -a); // v = a - b
	return v;
}

gsl_vector* I(gsl_vector* x)
{
	gsl_vector* v = gsl_vector_alloc(x->size);
	for (int i = 0; i < x->size; i++) {
		double val = gsl_vector_get(x, i);
		if (val > 0) {
			gsl_vector_set(v, i, val);
		}
		else {
			gsl_vector_set(v, i, 0);
		}
	}
	return v;
}

double E(gsl_vector* x, gsl_vector* p)
{
	/* double e = 0;
	for (int i = 0; i < x->size; i++) {
		e += (gsl_vector_get(x, i) * gsl_vector_get(p, i));
	}
	return e;
	*/
	return gsl_stats_wmean(p->data, 1, x->data, 1, x->size);
	
}

gsl_vector* p_norm(gsl_vector* x, double mu, double sigma)
{
	gsl_vector* v = gsl_vector_alloc(x->size);
	for (int i = 0; i < x->size; i++) {
		// use transformation x - mu, where x - mu ~ N(0, sigma)
		double val = gsl_vector_get(x, i) - mu;
		double p = 0;
		if (i == 0) {
			// first
			p = gsl_cdf_gaussian_P(val, sigma);
		}
		else if (i == x->size - 1) {
			// tail
			p = gsl_cdf_gaussian_P(val, sigma)
				- gsl_cdf_gaussian_P(val - 1, sigma)
				+ gsl_cdf_gaussian_Q(val, sigma);
		}
		else {
			// others
			p = gsl_cdf_gaussian_P(val, sigma)
				- gsl_cdf_gaussian_P(val - 1, sigma);
		}
		gsl_vector_set(v, i, p);
	}
	return v;
}


double relative_error(double m, double n)
{
	// m is the old model's value
	// n is the new model's value
	// the result shouldn't be negative
	return (n - m) / n;
}

gsl_vector* range(int begin, int end, int interval)
{
	int n = (end - begin) / interval + 1;
	if (interval == 1) {
		n--;
	}
	gsl_vector* v = gsl_vector_alloc(n);
	for (int x = begin, j = 0; x < end; x += interval, j++) {
		gsl_vector_set(v, j, x);
	}
	return v;
}
